function[V] = gradient_fun_y3_avg(x, Yt, vdjQt, t, retailer, side)

global D I xij_for_investor     

Yt      = sortrows(Yt, [1,16,17]); 
n       = size(Yt,1); 
idxE   = (Yt(:,2)==2);
NOTC =size(Yt(Yt(:,9)==1,:),1); 
shareOTC = NOTC./n;
mu      = x(1);
sigma  = x(2);
c         = x(3); 

idxj        = find(Yt(:,17)==1); 
dealer_ids  = Yt(idxj,16);
idxd        = find(ismember(D{t},dealer_ids));

sigma0   = Yt(1,15);
xij     = Yt(idxj,14);

if xij_for_investor==1
    xid     = 0; %set to 0 if psi does not include xid
    xij_shocks = Yt(:,14) ; %set to 0 if price excludes xij 
    
elseif  xij_for_investor==0
    xid     = xij; %set to 0 if psi does not include xid
    xij_shocks =  0; %Yt(:,14);  %set to 0 if price excludes xij 
end 
q       = Yt(idxj,18);

fun = @(x) normpdf(x,mu,sigma);
Fun = @(x) normcdf(x,mu,sigma);
funprime   = @(x) -(exp(-x.^2./2).*x)./(sqrt(2.*pi));
 
theta_long   = Yt(:,5); 
price        = Yt(:,8);  

if side==1
    shocks = price - theta_long + xij_shocks;
elseif side==-1
    shocks = price - theta_long - xij_shocks;
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INSTITUIONALS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if retailer==0

thetaj = zeros(I{t},1);
for i=1:I{t}
    thetaj(i)  = mean(Yt(Yt(:,16)==D{t}(i),5));
end 

psi =  NaN*ones(I{t},1);
if side==1
        psi(idxd,:) =  sigma0*log(sum(exp(1/sigma0*(xij+q))))-xid;
elseif side==-1
        psi(idxd,:) =  -sigma0*log(sum(exp(1/sigma0*(xij-q))))+xid;
end 

psi_long   = zeros(size(Yt,1),1);
theta_long = zeros(size(Yt,1),1);

for t=1:size(Yt,1)
    idx0            = find(Yt(t,16)==dealer_ids);
    psi_long(t)     = psi(idx0);
    theta_long(t)   = thetaj(idx0);
end 

if side==1
    ask_can=psi_long;
elseif side==-1
    bid_can=psi_long;
end 
theta = theta_long;

if side==1 
    
    b = Inf;
    a = (ask_can - theta - mu+c)./sigma;

    Mills = (fun(a))./(1 - Fun(a));

    ExpfunT      =  mu + sigma.*Mills; 
    VarfunT      =  sigma.^2 .*(1 + (a .*fun(a) ) ./(1 -  Fun(a)) - ( Mills).^2);
    
    massOTC  = 1 - Fun(a);
    massE    = Fun(a);
    probOTC  = massOTC./(massOTC+massE); 

    dExpfunMu    = 1 + (-fun((c - mu - theta + ask_can)./sigma).^2 + (-1 + Fun((c - mu - theta + ask_can)./sigma)) ...
                   .*funprime((c - mu - theta + ask_can)./sigma))./(-1 + Fun((c - mu - theta + ask_can)./sigma)).^2;
    
    dExpfunSigma = (fun((c - mu - theta + ask_can)./sigma).* (sigma - sigma .*Fun((c - mu - theta + ask_can)./sigma) ...
                   + fun((c - mu - theta + ask_can)./sigma).* (-c + mu + theta - ask_can)) + (-1 + Fun((c - mu - theta + ask_can)./sigma)) ...
                   .*funprime((c - mu - theta + ask_can)./sigma).* (c - mu - theta + ask_can))./(sigma .*(-1 + Fun((c - mu - theta + ask_can)./sigma)).^2);
    
    dExpfunC     = (fun((c - mu - theta + ask_can)./sigma).^2 - (-1 + Fun((c - mu - theta + ask_can)./sigma)) ...
                  .*funprime((c - mu - theta + ask_can)./sigma))./(-1 + Fun((c - mu - theta + ask_can)./sigma)).^2;
    
    dVarfunMu  =(2 .*sigma .*fun((c - mu - theta + ask_can)./sigma).^3 - (-1 + Fun((c - mu - theta + ask_can)./sigma)).^2 ...
                .*funprime((c - mu - theta + ask_can)./sigma).* (c - mu - theta + ask_can) + fun((c - mu - theta + ask_can)./sigma)...
                .*(-1 + Fun((c - mu - theta + ask_can)./sigma)).* (sigma - sigma .*(Fun((c - mu - theta + ask_can)./sigma) ...
                + 2 .*funprime((c - mu - theta + ask_can)./sigma)) + fun((c - mu - theta + ask_can)./sigma).* (c - mu - theta + ask_can)))./(1 - Fun((c - mu - theta + ask_can)./sigma)).^3;
            
    dVarfunSigma =  2.*sigma - (((-1 + Fun((c - mu - theta + ask_can)./sigma)).* (sigma.* fun((c - mu - theta + ask_can)./sigma) ...
                + funprime((c - mu - theta + ask_can)./sigma).* (-c + mu + theta - ask_can)) + fun((c - mu - theta + ask_can)./sigma).^2 ...
                .*(c - mu - theta + ask_can)) .*(2 .*sigma .*fun((c - mu - theta + ask_can)./sigma) + (-1 + Fun((c - mu - theta + ask_can)./sigma))...
                .*(c -  mu - theta + ask_can)))./(sigma .*(-1 + Fun((c - mu - theta + ask_can)./sigma)).^3);

    dVarfunC = (-2.* sigma .*fun((c - mu - theta + ask_can)./sigma).^3 + fun((c - mu - theta + ask_can)./sigma).* (-1 + Fun((c - mu - theta + ask_can)./sigma))...
               .*(sigma.* (-1 + Fun((c - mu - theta + ask_can)./sigma) + 2.* funprime((c - mu - theta + ask_can)./sigma)) + fun((c - mu - theta + ask_can)./sigma)...
               .*(-c + mu + theta - ask_can)) + (-1 + Fun((c - mu - theta + ask_can)./sigma)).^2 .*funprime((c - mu - theta + ask_can)./sigma) ...
               .*(c - mu - theta + ask_can))./(1 -Fun((c - mu - theta + ask_can)./sigma)).^3;
    
    dProbOTCcondMu = fun((c - mu - theta + ask_can)./sigma)./sigma;

    dProbOTCcondSigma = ((c - mu - theta + ask_can) .*fun((c - mu - theta + ask_can)./sigma))./sigma.^2;

    dProbOTCcondC = -fun((c - mu - theta + ask_can)./sigma)./sigma;
    

elseif side==-1 
     
    b = (bid_can -c - theta - mu)./sigma;
    a = -Inf;
    
    Mills = ( - fun(b))./(Fun(b) );

    ExpfunT      =  mu + sigma.*Mills; 
    VarfunT      =  sigma.^2 .*(1 + ( -  b.*fun(b)) ./(Fun(b) ) - ( Mills).^2);
    
    massOTC  = Fun(b) ;
    massE    = 1-Fun(b);
    probOTC  = massOTC./(massOTC+massE);
    
    dExpfunMu = 1 + (-fun(-((c + mu + theta - bid_can)./sigma)).^2 + Fun(-((c + mu + theta - bid_can)./sigma)) ...
                .*funprime(-((c + mu + theta - bid_can)./sigma)))./Fun(-((c + mu + theta - bid_can)./sigma)).^2;
    
    dExpfunSigma = (-sigma.*fun(-((c + mu + theta - bid_can)./sigma)).* Fun(-((c + mu + theta - bid_can)./sigma)) ...
                + fun(-((c + mu + theta - bid_can)./sigma)).^2 .*(c + mu + theta - bid_can) - Fun(-((c + mu + theta - bid_can)./sigma))...
                .*funprime(-((c + mu + theta - bid_can)./sigma)).* (c + mu + theta - bid_can))./(sigma.* Fun(-((c + mu + theta - bid_can)./sigma)).^2);
    
    dExpfunC = (-fun(-((c + mu + theta - bid_can)./sigma)).^2 + Fun(-((c + mu + theta - bid_can)./sigma)).* funprime(-((c + mu + theta - bid_can)./sigma)))...
             ./Fun(-((c + mu + theta - bid_can)./sigma)).^2;
    
     dVarfunMu = (1./(Fun(-((c + mu + theta - bid_can)./sigma)).^3)).*(-2.* sigma .*fun(-((c + mu + theta - bid_can)./sigma)).^3 ...
               + fun(-((c + mu + theta - bid_can)./sigma)).* Fun(-((c + mu + theta - bid_can)./sigma)).* (sigma.* (Fun(-((c + mu + theta - bid_can)./sigma)) ...
               + 2.* funprime(-((c + mu + theta - bid_can)./sigma))) + fun(-((c + mu + theta - bid_can)./sigma)).* (c + mu + theta - bid_can)) ...
               - Fun(-((c + mu + theta - bid_can)./sigma)).^2 .*funprime(-(( c + mu + theta - bid_can)./sigma)).* (c + mu + theta - bid_can));
   
    dVarfunSigma = 2 .*sigma + ((2 .*sigma .*fun(-((c + mu + theta - bid_can)./sigma)) - Fun(-((c + mu + theta - bid_can)./sigma))...
                .*(c + mu + theta - bid_can)).* (-sigma .*fun(-((c + mu + theta - bid_can)./sigma)).* Fun(-((c + mu + theta - bid_can)./sigma))...
                + fun(-((c + mu + theta - bid_can)./sigma)).^2 .* (c + mu + theta - bid_can) - Fun(-((c + mu + theta - bid_can)./sigma))...
                .*funprime(-((c + mu + theta - bid_can)./sigma)).* (c + mu + theta - bid_can)))./(sigma .*Fun(-((c + mu + theta - bid_can)./sigma)).^3);
    
     dVarfunC = (1./(Fun(-((c + mu + theta - bid_can)./sigma)).^3)).*(-2 .*sigma.* fun(-((c + mu + theta - bid_can)./sigma)).^3 ...
         + fun(-((c + mu + theta - bid_can)./sigma)) .*Fun(-((c + mu + theta - bid_can)./sigma)).* (sigma.* (Fun(-((c + mu + theta - bid_can)./sigma))...
         + 2 .*funprime(-((c + mu + theta - bid_can)./sigma))) + fun(-((c + mu + theta - bid_can)./sigma)).* (c + mu + theta -  bid_can))...
         - Fun(-((c + mu + theta - bid_can)./sigma)).^2 .*funprime(-((c + mu + theta - bid_can)./sigma)).* (c + mu + theta - bid_can));
        
     dProbOTCcondMu = -fun(-((c + mu + theta - bid_can)./sigma))./sigma;
     
     dProbOTCcondSigma= ((c + mu + theta - bid_can) .*fun(-((c + mu + theta - bid_can)./sigma)))./sigma.^2;
     
     dProbOTCcondC = -fun(-((c + mu + theta - bid_can)./sigma))./sigma;
     
end 
    
    moment1 = (ExpfunT - shocks);
    moment2 = (VarfunT -(shocks -ExpfunT).^2); 
    
    moment1(idxE)=0;
    moment2(idxE)=0;

    moment3 = probOTC - shareOTC;
    moment = [moment1, moment2, moment3];
    
    S = (moment'*moment);
    S=1/n*(S);

    
    dExpfunMu(idxE)=0;
    dExpfunSigma(idxE)=0;
    dExpfunC(idxE)=0;
    dVarfunMu(idxE)=0;
    dVarfunSigma(idxE)=0;
    dVarfunC(idxE)=0;

     G = [sum(dExpfunMu)/n,   sum(dVarfunMu)/n,    sum(dProbOTCcondMu)/n;...
        sum(dExpfunSigma)/n, sum(dVarfunSigma)/n, sum(dProbOTCcondSigma)/n;...
        sum(dExpfunC)/n,     sum(dVarfunC)/n,     sum(dProbOTCcondC)/n];
       
end    
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RETAILERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%

if retailer==1 
if side==1 
    
   Mills = 0;
   ExpfunT  =  mu + sigma.*Mills; 
   VarfunT  =  sigma.^2 ;

   dExpfunMu    = ones(n,1)*1 ;
   dExpfunSigma =ones(n,1)* 0;
   dVarfunMu    =ones(n,1)* 0;
   dVarfunSigma =ones(n,1)* 2 .*sigma;
       
elseif side==-1 

   Mills = 0;
   ExpfunT      =  mu + sigma.*Mills; 
   VarfunT      =  sigma.^2 ;
  
   dExpfunMu    = ones(n,1)*1;
   dExpfunSigma = ones(n,1)*0;
     
   dVarfunMu    = ones(n,1)*0;
   dVarfunSigma = ones(n,1)* 2 .*sigma;
   
end 

    moment1 = (ExpfunT - shocks);
    moment2 = (VarfunT -(shocks -ExpfunT).^2); 
    moment = [moment1, moment2];
    
    S = (moment'*moment);
    S=1/n*(S);

     G = [sum(dExpfunMu)/n,   sum(dVarfunMu)/n ;...
        sum(dExpfunSigma)/n, sum(dVarfunSigma)/n];
    
end 
    
V=inv(G*inv(S)*G');

end 